""" Crystal strucutres of the VASP's POSCAR-format files. """
import sys        # call sys.exit() at errors.
import os         # call os.getcwd(); os.sep.
import os.path    # os.path.isfile()
import numpy as np

# A helper function that converts True to 'T' and False to 'F'
def bool_to_char_(t):
    if t : return 'T'
    else : return 'F'

# Chemical elements
ELEMENTS = {
        1: ("H",  "hydrogen"),
        2: ("He", "helium"),
        3: ("Li", "lithium"),
        4: ("Be", "beryllium"),
        5: ("B",  "boron"),
        6: ("C",  "carbon"),
        7: ("N",  "nitrogen"),
        8: ("O",  "oxygen"),
        9: ("F",  "fluorine"),
        10: ("Ne", "neon"),
        11: ("Na", "sodium"),
        12: ("Mg", "magnesium"),
        13: ("Al", "alluminum"),
        14: ("Si", "silicon"),
        15: ("P",  "phosphorous"),
        16: ("S",  "sulphur"),
        17: ("Cl", "chlorine"),
        18: ("Ar", "argon"),
        19: ("K",  "potassium"),
        20: ("Ca", "calsium"),
        21: ("Sc", "scandium"),
        22: ("Ti", "titanium"),
        23: ("V",  "vanadium"),
        24: ("Cr", "chromium"),
        25: ("Mn", "manganese"),
        26: ("Fe", "iron"),
        27: ("Co", "cobalt"),
        28: ("Ni", "nickel"),
        29: ("Cu", "copper"),
        30: ("Zn", "zinc"),
        31: ("Ga", "gallium"),
        32: ("Ge", "gemanium"),
        33: ("As", "arsenic"),
        34: ("Se", "selenium"),
        35: ("Br", "bromine"),
        36: ("Kr", "krypton"),
        37: ("Rb", "rubidium"),
        38: ("Sr", "strontium"),
        39: ("Y",  "yttrium"),
        40: ("Zr", "zirconium"),
        41: ("Nb", "niobium"),
        42: ("Mo", "molybdenum"),
        43: ("Tc", "technetium"),
        44: ("Ru", "rhodium"),
        45: ("Pd", "paladium"),
        46: ("Ag", "silver"),
        47: ("Cd", "cadmium"),
        48: ("In", "indium"),
        49: ("Sn", "tin"),
        50: ("Sb", "antimony"),
        52: ("Te", "tellurium"),
        53: ("I",  "iodine"),
        54: ("Xe", "xenon"),
        55: ("Cs", "cesium"),
        56: ("Ba", "barrium"),
        57: ("La", "lanthanum"),
        58: ("Ce", "cerium"),
        59: ("Pr", "praseodymium"),
        60: ("Nd", "neodymium"),
        61: ("Pm", "promethium"),
        62: ("Sm", "samarium"),
        63: ("Eu", "europium"),
        64: ("Gd", "gadolinium"),
        65: ("Tb", "terbium"),
        66: ("Dy", "dysprosium"),
        67: ("Ho", "holmium"),
        68: ("Er", "erbium"),
        69: ("Tm", "thulium"),
        70: ("Yb", "ytterbium"),
        71: ("Lu", "lutetium"),
        72: ("Hf", "hafnium"),
        73: ("Ta", "tantalum"),
        74: ("W",  "tungsten"),
        75: ("Re", "rhenium"),
        76: ("Os", "osmium"),
        77: ("Ir", "iridium"),
        78: ("Pt", "platinum"),
        79: ("Au", "gold"),
        80: ("Hg", "mercury"),
        81: ("Tl", "thallium"),
        82: ("Pb", "lead"),
        83: ("Bi", "bismuth"),
        84: ("Po", "polium"),
        85: ("At", "astatine"),
        86: ("Rn", "radon"),
        }

class VASP_structure:
    """ VASP format strcutres.
    
    A structure contains comment, scalar, basis vectors,
    format (vasp4 or vasp5), selective dynamics flag ant tags, 
    names of species, numbers of atoms of the species, and atom positions.
    """
    def __init__(self):
        """ Initializes an empty, but valid, structure."""
        self.comment = None
        self.scalar = 1.0
        self.bv = np.array([[1.0, 0, 0], [0.0, 1.0, 0.0], [0.0, 0, 1.0]])
        self.species = []
        self.n_atoms = []
        self.format = None      # VASP4 or VASP5
        self.sd_flag = None     # selective dynamics
        self.sd_tags = []       # an array such as [[True, True, True], ...]
        self.atom_positions = []

        # Inversion of bv. This will be updated whenever bv is modified
        self.bvi_ = np.array([[1.0, 0, 0], [0.0, 1.0, 0.0], [0.0, 0, 1.0]])

    def read(self, POSCAR):
        """initialize a structure from a POSCAR-style file.

        POSCAR-style files include POSCAR, CONTCAR, CHG, CHGCAR, XDATCAR.
        NOTE: XDATCAR contains multiple structures. This method is to read
        the first structure. Another method is to be defined for reading
        all structures.
        
        Arguments:
        POSCAR: a file name or a ready-to-read input stream.
            If POSCAR is a file name, a file will be opened read-only, 
            then read, and closed. If POSCAR is an opened file,
            it will not be closed in this function; the caller is
            responsible to close the file.
        Returns:
            void. The structure is updated.
        Exceptions:
            IOError, ValueError.
        """
        # if POSCAR is a string, open a file read-only
        if type(POSCAR) == str:
            if POSCAR[0]=='/':  # POSCAR is an absolute path.
                F_POSCAR = open(POSCAR, "r")
            else:   # POSCAR is a relative path.
                F_POSCAR = open(os.getcwd()+os.sep+POSCAR, "r")
        else: 
            F_POSCAR = POSCAR  # POSCAR is a file that is already opend.

        self.comment = F_POSCAR.readline().strip()
        self.scalar = float(F_POSCAR.readline())
        # bv is columnwise: each column is a basis vector
        self.bv[:, 0] = [float(x) for x in F_POSCAR.readline().split()]
        self.bv[:, 1] = [float(x) for x in F_POSCAR.readline().split()]
        self.bv[:, 2] = [float(x) for x in F_POSCAR.readline().split()]
        self.bv *= self.scalar

        # Because bv is updated, now updates bvi_.
        # If bv is not valid (e.g., some vectors are coplanar, etc.), this
        # operation will raise an error.
        self.bvi_ = np.array((np.matrix(self.bv)*self.scalar).I)

        # Determines whether it is VASP5 or VASP4 format.
        element_or_number = F_POSCAR.readline().strip()
        if(element_or_number[0].isalpha()):
              self.format = "VASP5"
              self.species = element_or_number.split()
              self.n_atoms = [int(x) for x in F_POSCAR.readline().split()]
              if len(self.species) != len(self.n_atoms):
                   raise ValueError(
                       "numbers of species and atom-types are unequal!"
                   )
        else:
             self.format = "VASP4"
             self.n_atoms = [int(x) for x in element_or_number.split()]
             for i, x in enumerate(self.n_atoms):
                 self.species.append("UN"+chr(i+65))
                 # chr(65) is A, chr(66) is B, and so on.
                 # TODO: when there are POTCAR, OUTCAR, or vasprun.xml,
                 # it is possible to deduce the names of the elements.

        # Determines whether selective-dynamics is used.
        LINEBUFFER = F_POSCAR.readline()
        if LINEBUFFER[0] in "sS":
            self.sd_flag = True
            LINEBUFFER = F_POSCAR.readline()
        else:
            self.sd_flag = False
        
        # Determines which coordinate system is used.
        coordinate_system = -1    # 0: direct; 1: Cartesian.
        if LINEBUFFER[0] in "cCkK":
            coordinate_system = 1   # Cartesian coordinate system
        else:
            coordinate_system = 0   # Direct (or fractional) coordinate system

        # Reads in different fields, depending on selective-dynamics flag
        n_total = np.sum(self.n_atoms)
        if self.sd_flag == False:
            if coordinate_system==0:
                for i in range(n_total):
                    x, y, z = [float(t) for t in F_POSCAR.readline().split()]
                    # Corrects out-of-box coordinates
                    x -= np.floor(x)
                    y -= np.floor(y)
                    z -= np.floor(z)
                    self.atom_positions.append([x, y, z])
                    self.sd_tags.append([False, False, False])  # not used.
            elif coordinate_system==1:
                # Always store coordinates internally in fractional format.
                for i in range(n_total):
                    c = [float(t) for t in F_POSCAR.readline().split()]
                    self.atom_positions.append(
                        internal_coordinate_(self.bvi_, c)
                    )
                    self.sd_tags.append([False, False, False])  # not used.
            else: 
                raise ValueEror("Unknown coordinate system.")
        elif self.sd_flag == True:
            if coordinate_system==0:
                for i in range(n_total):
                    m = F_POSCAR.readline().split()
                    x, y, z = [float(t) for t in m[0:3] ]
                    # Corrects out-of-box coordinates
                    x -= np.floor(x)
                    y -= np.floor(y)
                    z -= np.floor(z)
                    self.atom_positions.append([x, y, z])

                    # This syntax is not allowed in python:
                    #t1, t2, t3 = [t[0] in 'tT' ? True : False for t in m[3:6]]
                    if m[3][0] in "tT": t1 = True
                    else: t1 = False
                    if m[4][0] in "tT": t2 = True
                    else: t2 = False
                    if m[5][0] in "tT": t3 = True
                    else: t3 = False
                    self.sd_tags.append([t1, t2, t3])
            elif coordinate_system==1:
                # Always store coordinates internally in fractional format.
                for i in range(n_total):
                    m = F_POSCAR.readline().split()
                    c = [float(t) for t in m[0:3] ]
                    self.atom_positions.append(
                        internal_coordinate_(self.bvi_, c)
                    )

                    if m[3][0] in "tT": t1 = True
                    else: t1 = False
                    if m[4][0] in "tT": t2 = True
                    else: t2 = False
                    if m[5][0] in "tT": t3 = True
                    else: t3 = False
                    self.sd_tags.append([t1, t2, t3])
            else:
                raise ValueError("Unknown coordinate system.")
        else:
            raise ValueError("Unknown tags for selective dynamics.")

	# Closes file POSCAR if it was opened in this function.
        if type(POSCAR) == str:
                F_POSCAR.close()


    def atom_name(self, i):
        """ Returns the element name of index i"""
        n_total = np.sum(self.n_atoms)
        if i>=n_total:
            print("VASP_structure.atom_name(): index exceeds total number.")
            raise ValueError
        if i<0:
            print("VASP_structure.atom_name(): index is incorrect.")
            raise ValueError
        accumulate_n1 = 0
        accumulate_n2 = 0
        for ind, n in enumerate(self.n_atoms):
            accumulate_n2 += n
            if i<accumulate_n2 and i>=accumulate_n1:
                return self.species(ind)
            else:
                accumulate_n1 = accumulate_n2


    def addAtomAtFractionalCoordinate(self, species, position):
        """ add an atom at a given fractional coordinate.

        Arguments:
            species: (string) name of the atom type, such as "Co", "Fe1", etc.
                If the species is already present in the structure, then this
                atom is appended to the list of the existent species; 
                otherwise a new list is created to contain this atom.
            position: (3-member list) fractional coordinate of the newly added
                atom, such as [0.25, 0.25, 0.75]. It can also be a numpy
                array. If any component is outsiede the range [0,1), it will
                be displaced according to the periodic boundary condition
                into the range [0, 1).

        Returns:
            void. The structure is updated.

        Note about the species:
        if the format is VASP4, the species are "UNA", "UNB", "UNC", etc.;
        if the format is VASP5, the species are abbreviation of element names.

        TODO: the new atom is always appended at the end of the species.
        It may require the flexibility to insert the atom at a given index.
        """
        # x,y,z should be normalized into the range [0, 1)
        c_internal = [x-np.floor(x) for x in position]

        if species in self.species:
            # Gets the index of the species
            i = self.species.index(species)
            # Counts the number of atoms of species before this index
            n_pre = np.sum(self.n_atoms[0:i])
            n_cur = self.n_atoms[i]
            self.atom_positions.insert(n_pre+n_cur, c_internal)
            self.sd_tags.insert(n_pre+n_cur, [False, False, False])
            self.n_atoms[i] = n_cur + 1
        else:
            # appends a specie, and a number of species
            self.species.append(species)
            self.n_atoms.append(1)
            self.atom_positions.append(c_internal)
            self.sd_tags.append(n_pre+n_cur, [False, False, False])


    def addAtomAtCartesianCoordinate(self, species, position, index=-1):
        """ add an atom at a given Cartesian coordinate.

        Arguments:
            species: (string) name of the atom type, such as "Co", "Fe1", etc.
                If the species is already present in the structure, then this
                atom is appended to the list of the existent species; 
                otherwise a new list is created to contain this atom.
            position: (3-number list) Cartesion coordinate of the newly added
                atom, such as [10.01, 7.26, -5.32]. This coordinate will be
                be transformed into internal fractional coordinate.

        Returns:
            void. The structure is updated.

        Note about the species:
        if the format is VASP4, the species are "UNA", "UNB", "UNC", etc.;
        if the format is VASP5, the species are abbreviation of element names.
        """
        c_internal = internal_coordinate_(self.bvi_, position)
        self.addAtomAtFractionalCoordinate(species, c_internal)


    def expandCell(self, s, ref_center=(0,0,0)):
        """ changes lattice constants (a,b,c), keeps atom positions unchanged.

        Argument:
            s: (positive float scalar) ratio between new lattice constant
                and old lattice constant. When s<1.0, the actual result is 
                to shrink the cell. Examples are 1.2, 25.0/20.0, 0.8.
            ref_center: (3-number list) the fractional coordinate around which
                the cell is expanded (or shrinked). This cooridnate will be 
                unchanged before and after the transformation.

        Returns:
            void. The structure is updated.
        """
        # s may be a scalar or a 3-member list.
        if type(s) is list:
            if len(s)!=1 or len(s)!=3:
                print("In expandCell(): scale should be 1- or 3-number array.")
                raise ValueError
            if any(x<0.0 for x in s):
                print("In expandCell(): scale factor should be positive.")
                raise ValueError
            if len(s)==1:
                self.bv *= s[0]
                sr_x = 1.0/s[0]
                sr_y = 1.0/s[0]
                sr_z = 1.0/s[0]
            elif len(s)==3:
                self.bv = np.matmul(self.bv, np.diag(s))
                sr_x = 1.0/s[0]
                sr_y = 1.0/s[1]
                sr_z = 1.0/s[2]
            else:
                pass
        else: # s need to be a positive number
            if s<=0.0:
                print("In expandCell(): scale factor", s, "should be positive.")
                print("This is unreasonable. I will stop here.")
                raise ValueError
            self.bv *= s
            sr_x = 1.0/s
            sr_y = 1.0/s
            sr_z = 1.0/s

        for i, c in enumerate(self.atom_positions):
            temp_x = (c[0]-ref_center[0]) * sr_x + ref_center[0]; 
            temp_y = (c[1]-ref_center[1]) * sr_y + ref_center[1]; 
            temp_z = (c[2]-ref_center[2]) * sr_z + ref_center[2]; 
            self.atom_positions[i][0] = temp_x
            self.atom_positions[i][1] = temp_y
            self.atom_positions[i][2] = temp_z



    def set_sd(self, sd_flag=True):
        """ Applies selective dynamic flag and tags.

        if the structure already has sd_flag=True, this function will show
        a warning message, and does nothing.

        Arguments:
            sd_flag: (boolean) Whether to apply the selective dynamic tag.
                If this argument is True, and the structure already has 
                self.sd_flag=True, this call does nothing.
                If this argument is True, and the structure has 
                self.sd_flag=False, this call initialize a all-True sd_tags.
                If this argument is false, the structure's sd_tags will be
                reset to an all-True array.
                
        Returns:
            void. The structure is updated.
        """
        if sd_flag:
            if self.sd_flag :
                return
            else :
                self.sd_flag = True
                self.sd_tags = np.array(np.ones((np.sum(self.n_atoms), 3), 
                    dtype=np.int), dtype=np.bool).tolist()
        else :
            self.sd_tags = np.array(np.ones((np.sum(self.n_atoms), 3),
                    dtype=np.int), dtype=bool).tolist()


    def set_sd_tag(self, index , sd_tag):
        """ Sets a selective dynamic tag to an atom of the specified index.

        Arguments:
            index: (integer) the index of the atom to be applied the tag.
            sd_tag: (3-member list of boolean type) the selective dyanmic tag.
                Example: [True, False, False], [1, 0, 0].

        Returns:
            void. The structure is updated.
        """
        converted_bool_array = np.array(sd_tag, dtype=bool)
        self.sd_tags[index][0] = converted_bool_array[0]
        self.sd_tags[index][1] = converted_bool_array[1]
        self.sd_tags[index][2] = converted_bool_array[2]


    def write(self, POSCAR, selective_dynamics=True, coordinate_type='d'):
        """Writes the structure to a file.

        Arguments:
            POSCAR: (string) The output file's name, such as "POSCAR-new".
                If the file exists, the program asks the user whether to 
                overwrite.
            selective_dynamics: (boolean) If true, and if sd_flag of the 
                structure is also true, write tags of selective dynamics. 
                Otherwise (selectrive_dynamics is false, or sd_flag is false),
                no tags of selective dynamics are written.
            coordinate_type: (string) According to VASP's rule, if the first
                character is 'd' or 'D' (means 'Direct'), then coordinate type 
                is fractional; if the first character is 'c' or 'C' (menas
                cartesian), then coordinate type is cartesian. Other
                settings throw a ValueError exception.
        Returns:
            void. The output file is filled with the structure.
        """
        # Check the coordinate type
        if coordinate_type.casefold()[0]=='d' :
            CT = 0
        elif coordinate_type.casefold()[0]=='c' :
            CT = 1
        else :
            print("***Error: unknown coorindate type:", coordinate_type)
            raise ValueError

        # Check if the file exists; warn if the file exists.
        if os.path.isfile(POSCAR):
            print("***Warning: File " + POSCAR + " already exists.\n",
                  "  Action required. Overwrite it or not?\n"
                  "  Answer Y or y to overwrite, others to stop.\n")
            a = raw_input("Your answer: ")
            if a=="y" or a=="Y":
                F_POSCAR = open(POSCAR, "w")
            else:
                print("   Notify: The progrom stops, and no files have been modified.")
                return 1
        else:
            F_POSCAR = open(POSCAR, "w")

        # writes the content of the structure
        F_POSCAR.write(self.comment + "\n")
        F_POSCAR.write(str(self.scalar) + "\n")
        for i in range(3):
            F_POSCAR.write(" ".join([str(x) for x in self.bv[i,:]]) + "\n")
        if self.format == "VASP5" :
            F_POSCAR.write(" ".join(self.species) + "\n")
        F_POSCAR.write(" ".join([str(x) for x in self.n_atoms]) + "\n")

        # Writes positions and selective dynamic tags
        if selective_dynamics and self.sd_flag:
            F_POSCAR.write("Selective\n")
            if CT==0 :  # coordinate type is direct.
                F_POSCAR.write("Direct\n")
                for i, a in enumerate(self.atom_positions):
                    F_POSCAR.write(" ".join([str(x) for x in a]) + "   " )
                    F_POSCAR.write(" ".join(map(bool_to_char_, self.sd_tags[i]))
                        + "\n")
            else :    # coordinate type is cartesian
                F_POSCAR.write("Cartesian\n")
                for i, a in enumerate(self.atom_positions):
                    b = self.cartesian_coordinate(a)
                    F_POSCAR.write(" ".join([str(x) for x in b]) + "   " )
                    F_POSCAR.write(" ".join(map(bool_to_char_, self.sd_tags[i]))
                        + "\n")
        else:
            if CT==0 :  # coordinate type is direct.
                F_POSCAR.write("Direct\n")
                for a in self.atom_positions:
                    F_POSCAR.write(" ".join([str(x) for x in a]) +"\n" )
            else :    # coordinate type is cartesian
                F_POSCAR.write("Cartesian\n")
                for i, a in enumerate(self.atom_positions):
                    b = self.cartesian_coordinate(a)
                    F_POSCAR.write(" ".join([str(x) for x in b]) + "   " )
        F_POSCAR.close()


    def distance(self, index1, index2):
        """ Calculates distance between atoms index1 and index2.
        
        Arguments:
            index1, index2: (positive integer) indices of two atoms.
        
        Returns:
            The distance (float) between the two atoms. Periodic
                boundary conditions are considered, so that for example
                the difference between 0.01 and 0.99 is not 0.98, but 0.02.
        """
        # check whether index1 and index2 are valid
        n_total = np.sum(self.n_atoms)
        if index1<0 or index1>=n_total:
            print("***Error: atom index", index1, "is out of range")
            sys.exit(2)
        if index2<0 or index2>=n_total:
            print("***Error: atom index", index2, "is out of range")
            sys.exit(2)
        
        # get the fractional coordinates
        c1 = self.atom_positions[index1]
        c2 = self.atom_positions[index2]

        # apply the periodic boundary condition
        d = np.array(c1)-np.array(c2)
        # Well-structured model should have each component of d in [-1, 1]
        if d[0]<-0.5: d[0] += 1.0
        elif d[0]>0.5: d[0] -= 1.0
        else: pass
        if d[1]<-0.5: d[1] += 1.0
        elif d[1]>0.5: d[1] -= 1.0
        else: pass
        if d[2]<-0.5: d[2] += 1.0
        elif d[2]>0.5: d[2] -= 1.0
        else: pass

        # calculate the distance in Cartesian coordinate system.
        return magnitude_(self.bv, d)


    def cartesian_coordinate(self, r):
        """ Returns the cartesian coordinate of a fractional coordinate.
        
        r is an internal coordinate in the coordinate system of the structure.
        
        Arguments:
        r: (3-number ist or a numpy.ndarray) A Cartesian coordinate.
        
        Returns:
            The cartesian coordinate (a 3-number list).
        """
        return coordinate_in_system(self.bv, r)


    def NN(self):
        """ Returns nearest neighbors of every atom.

        Returns:
            a tuple of (atom_index, distance).
                The atom_index is a two-dimensial list of integer, n_total by n_total-1.
                The first index is for atoms, the second index is for neibours.
                The distance is also a two-dimensional list, n_total by n_total-1.
                The first index is for atoms, the second index is for distances.
        """
        n_total = np.sum(self.n_atoms)
        nn_index = []
        nn_distance = []
        for i in range(n_total): 
            nn_index.append([])
            nn_distance.append([])

        for i, r1 in enumerate(self.atom_positions) :
            for j, r2 in enumerate(self.atom_positions):
                if (i==j) or (i in nn_index[j]) or (j in nn_index[i]): 
                    continue
                d = self.distance(i, j)
                nn_index[i].append(j)
                nn_distance[i].append(d)
                nn_index[j].append(i)
                nn_distance[j].append(d)
        # sort according distance, and store in new arrays
        nn_index2 = []
        nn_distance2 = []
        for i, d in enumerate(nn_distance):
            nn_index2.append([])
            nn_distance2.append([])
            p = np.argsort(d)
            for m in p:
                nn_index2[i].append(nn_index[i][m])
                nn_distance2[i].append(nn_distance[i][m])

        return (nn_index2, nn_distance2)


    def corner2center(self):
        """ Changes the center of mass from a corner to cell's center.

        This is useful for molecular structures whose mass center is (0,0,0).
        """
        for p in self.atom_positions:
            if p[0]>=0.5: p[0] -= 0.5
            else: p[0] += 0.5
            if p[1]>=0.5: p[1] -= 0.5
            else: p[1] += 0.5
            if p[2]>=0.5: p[2] -= 0.5
            else: p[2] += 0.5
        

    def atomCenter(self):
        """ Calculates the center of atom positions (not mass center).

        This center does not necessarily coicide with the mass center,
        mainly because of the periodic boundary condition.

        Returns:
        a list [x, y, z]: Center of atom positions, in cartesian coordinate.
        """
        n_total = np.sum(self.n_atoms)
        center_x = np.sum([l[0] for l in self.atom_positions])/n_total
        center_y = np.sum([l[1] for l in self.atom_positions])/n_total
        center_z = np.sum([l[2] for l in self.atom_positions])/n_total
        center = self.cartesian_coordinate([center_x, center_y, center_z])
        return center
        

    def polygonCenter(self, p, n):
        """ Returns center of a polygon.
        
        A polygon is defined as a list of indices. Polygon p has n vertices. 
        This function is used to define pentagonCenter() and hexagonCenter().
    
        Arguments:
        p: [[x0, y0, z0], [x1, y1, z1], ..., [x{n-1}, y{n-1}, z{n-1}]].
        n: number of vertices of the polygon.
    
        Returns:
        [center_x, center_y, center_z]
        """
        x=0.0; y=0.0; z=0.0;
        for i in p:
            t = self.atom_positions[i]
            x += t[0]
            y += t[1]
            z += t[2]
        coordinate_fractional = (x/n, y/n, z/n)
        return self.cartesian_coordinate(coordinate_fractional)

    def pentagonCenter(self, p):
        """ Returns center of a pentagon. See polygeonCenter() for help."""
        return self.polygonCenter(p, 5)

    def hexagonCenter(self, p):
        """ Returns center of a hexagon. See polygeonCenter() for help."""
        return self.polygonCenter(p, 6)


    def stat(self):
        """ Prints statistics of the structure. For debugging."""
        print("Comment:", self.comment)
        print("Universal scalar:", self.scalar)
        print("Basis vector:")
        for i in range(3):
            print("\t{0:9.7f} {1:9.7f} {2:9.7f}".format(
                  self.bv[i][0], self.bv[i][1], self.bv[i][2]))
        print("Atom species:", self.species)
        print("Number of atoms of each species:", self.n_atoms)
        print("Total number of atoms:", np.sum(self.n_atoms))
        print("Structural format:", self.format)
        if self.sd_flag:
            print("Selective dynamics: Yes")
        else:
            print("Selective dynamics: No")


# Helper functions
def coordinate_in_system(b, v):
    """ Returns coordinate in a coordinate system."""
    x = b[0,0]*v[0] + b[0,1]*v[1] + b[0,2]*v[2]
    y = b[1,0]*v[0] + b[1,1]*v[1] + b[1,2]*v[2]
    z = b[2,0]*v[0] + b[2,1]*v[1] + b[2,2]*v[2]
    return (x, y, z)


def internal_coordinate_(bi, v):
    """ Calculate fractional coordinate of an absolute coordinate."""
    x, y, z = coordinate_in_system(bi, v)
    x -= np.floor(x)
    y -= np.floor(y)
    z -= np.floor(z)
    return [x, y, z]


def magnitude_(b, v):
    """ Calculate length of a vector in fractional coordinate"""
    x, y, z = coordinate_in_system(b, v)
    return np.sqrt(x*x + y*y + z*z)

def distances(s1, s2):
    """ returns the interatomic distances between two structures.

    The two structures s1 and s2 should be similar in that (a) their
    basis vectors are the same, and (b) numbers of atoms are the same.
    Otherwise, calculating interatomic distances are not meaningful.
    """
    diff_array = []
    for i, a1 in s1.atom_positions :
        a2 = s2.atom_positions[i]
        dx = a1[0] - a2[0]
        dy = a1[1] - a2[1]
        dz = a1[2] - a2[2]
        dist = np.sqrt(dx*dx + dy*dy + dz*dz)
        diff_array.append([dist, dx, dy, dz])
    return diff_array

# test code
if __name__=="__main__":
    s = VASP_structure()
    f_CONTCAR =  open("CONTCAR", "r")
    s.read(f_CONTCAR)
    # Show statitistics of the structure
    s.stat()

    f_CONTCAR.close()

